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ABSTRACT 


Galactic outflows driven by star formation or active galactic nuclei are typically formed by multi-phase gas whose temperature spans 
over 4 orders of magnitude. Probing the different outflow components requires multi-wavelength observations and long exposure 
times, especially in the distant Universe. So far, most of the high-z studies have focused on a single gas phase, but this kind of analysis 
may potentially miss a non-negligible fraction of the total outflowing gas content. In this work, we analyze the spatially resolved rest- 
frame UV and optical emission from HZ4, the highest redshift main sequence star-forming galaxy having a detected [Cn] outflow, 
which traces the neutral gas component. Our goal is to study the ionized interstellar medium in the galaxy and the properties of 
the ionized outflow as traced by the [Om]45007À and Ha emission lines. We exploit JWST/NIRSpec observations in the integral 
field spectroscopy mode to investigate the galaxy properties by making use of the brightest rest-frame optical emission lines. Their 
high spectral and spatial resolution allows us to trace the ionized outflow from broad line wings and spatially resolve it. We also 
re-analyze the [Ci] ALMA data to compare the neutral atomic and ionized outflow morphologies, masses, and energetics. We find 
that the system consists of a galaxy merger, instead of a rotating disk as originally inferred from low-resolution [Cn] observations, 
and hosts an extended ionized outflow. The ionized outflow is being launched from a region hosting an intense burst of star formation 
and extends over 4 kpc from the launch site. The neutral and ionized outflows are almost co-spatial, but the mass loading factor (the 
ratio between the mass outflow rate and the galaxy star formation rate) in the ionized gas phase is two orders of magnitude smaller 


1. Introduction 


Several theoretical models suggest that galactic outflows driven 
by star formation (SF) and active galactic nuclei (AGN) are cru- 
cial to explain the lack of galaxies in both the high- and low- 
mass ends of the galaxy mass function (e.g., Fabian 2012; Naab 
& Ostriker 2017; Somerville & Davé 2015), the inefficiency 
of galaxies in turning baryons into stars (Dekel & Silk 1986; 
Behroozi et al. 2019), the metal enrichment of the circumgalac- 
tic medium (CGM) and the intergalactic medium (IGM) (Oppen- 
heimer & Davé 2006), and the shape of the mass-metallicity re- 
lation (Chisholm et al. 2018). In particular, low-mass (« 10 Mo) 
galaxies are believed to self-regulate the buildup of their stellar 
mass through the action of SF-driven outflows (Hopkins et al. 
2012). The energy and momentum injected by young massive 
stars in the shape of stellar winds, supernova (SN) explosions, 
and radiation pressure can provide feedback, by heating and ex- 
pelling the gas from the galaxy (e.g., Somerville & Davé 2015). 
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than in the neutral phase, as found for other lower redshift multi-phase outflows. 


Key words. Galaxies: high-redshift — ISM: jets and outflows — ISM: kinematics and dynamics 


Local observations of low mass galaxies have revealed that 
they also host accreting supermassive black holes at their cen- 
ters (e.g., Sartori et al. 2015; Reines et al. 2020), indicating that 
SF activity can be regulated by AGN-driven outflows. Theoret- 
ical works and simulations have indeed pointed out that AGN 
can boost the outflow temperature and velocity by two orders 
of magnitude (Silk 2017; Dashyan et al. 2018; Koudmani et al. 
2019), which means that AGN-driven outflows may affect the 
CGM around low-mass galaxies by polluting it with metals and 
heating it. Koudmani et al. (2021) have also shown that AGN- 
feedback in low-mass galaxies may be more relevant at z > 2, 
when AGN are expected to be more active. In conclusion, the 
primary mechanisms responsible for regulating star formation 
in low-mass galaxies are still unclear and debated. Theoretical 
models offer a wide range of solutions to this problem, and multi- 
wavelength observations are needed to test them. 


In the last ten years, observations and simulations have 
shown that galactic outflows include gas at different tempera- 
tures, densities, and states (molecular, neutral atomic, and ion- 
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ized) (e.g., Muratov et al. 2015; Janssen et al. 2016; Li et al. 
2017; Nelson et al. 2019; Fluetsch et al. 2019, 2021). In par- 
ticular, four components have been identified: a) the hot (T 
~ 10677 K) ionized outflowing gas, which emits in the X-rays 
band; b) the warm ionized component (T ~ 10^ K), which is usu- 
ally detected as blueshifted broad wings in the brightest optical 
emission lines, such as [Om]45007À and Ho; c) the cool neu- 
tral atomic part (T ~ 100 K) that is usually observed via NaID 
445890,5896À in the local Universe and [Ci]A158um! in high-z 
galaxies; and d) the cold molecular component (T~ 10 K), traced 
usually by the broad wings of the CO molecular transitions or 
by P-Cygni profiles of molecular transitions (e..g OH 119 um, 
H20). Aside from a few studies, the majority of the outflows are 
observed only through one tracer, making it difficult to estimate 
the global outflow properties such as the total outflow mass and 
mass loading factor, hence the effect on the host galaxy. 

While at cosmic noon (1 x z < 3) outflows are usually iden- 
tified from broad, often blueshifted components of optical emis- 
sion lines (e.g., Coatman et al. 2019; Förster Schreiber et al. 
2019; Villar Martín et al. 2020; Concas et al. 2022; Cresci et al. 
2023), at high redshift (z > 3) outflow studies have focused 
on the far infrared (FIR) emission lines. This is because optical 
emission lines are redshifted to near-infrared (NIR) wavelengths, 
which become hard to access from the ground. In particular At- 
acama Large Millimetre Array (ALMA) observations have re- 
vealed outflow signatures in molecular transitions (P-Cygni pro- 
files of OH 119 um, broad CO and H50 line wings; Herrera- 
Camus et al. 2019; Jones et al. 2019; Butler et al. 2023) and 
the atomic carbon transition, [Cr] (Gallerani et al. 2018; Ginolfi 
et al. 2020; Bischetti et al. 2019; Herrera-Camus et al. 2021; 
Solimano et al. 2024; Tripodi et al. 2023). However, the major- 
ity of these studies have focused mainly on AGN host galaxies, 
massive sub-millimeter galaxies, and lensed dusty star-forming 
galaxies (Jones et al. 2019; Spilker et al. 2020). Evidence for 
SF-driven neutral outflows in galaxies in the early Universe has 
been observed for the first time by stacking the [Cri] emission of 
samples of galaxies at z ~ 4 — 5 (Gallerani et al. 2018; Ginolfi 
et al. 2020). So far, broad [Cu] in the spectra of individual main- 
sequence galaxies has been reported only for two z > 4 sources 
(Herrera-Camus et al. 2021; Solimano et al. 2024). 

With the advent of JWST, we finally have access to the rest- 
frame optical emission lines of high-redshift galaxies, which al- 
lows us to study the outflows of high-z galaxies with the same 
tools developed for low-z ones. This also gives the possibility to 
simultaneously characterise up to two (or three, see above) out- 
flow components in the same target. The first studies with JWST 
have revealed that a large fraction of star-forming galaxies have 
ionized outflows, up to the cosmic dawn (Carniani et al. 2023; 
Zhang et al. 2024; Calabró et al. 2024). Therefore, the synergy 
between ALMA and JWST promises to be a rich probe into the 
physical nature of high-z outflows. 

In this work, we focus on HZ4 (RA = 09h58m28.5s, Dec= 
+02d03m06.7s), which is the most distant star-forming galaxy 
with evidence of neutral outflows through the [Cn] emission 
line (Herrera-Camus et al. 2021). HZ4 was identified as a Ly- 
man Break Galaxy (LBG) at z ~ 5.5 in the Cosmic Evolution 
Survey (COSMOS) field (Scoville et al. 2007) and was spec- 
troscopically confirmed with Keck Deep Extragalactic Imaging 
Multi-Object Spectrograph (DEIMOS) observations (Mallery 
et al. 2012) which detected Lya emission (Cassata et al. 2020; 
Le Févre et al. 2020). It was then followed up with ALMA 


! The fraction of [Cri] emission arising from ionized gas is only 1-30% 
of the total emission (Díaz-Santos et al. 2017; Cicone et al. 2018) 
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observations to target the [Cu] and dust continuum emission 
with low- (1.01"x 0.85") and high- (0.39"x 0.34") angular 
resolution (Capak et al. 2015; Béthermin et al. 2020; Herrera- 
Camus et al. 2021; Jones et al. 2021). In the ALPINE survey 
(Le Févre et al. 2020; Béthermin et al. 2020), HZ4 1s dubbed 
DEIMOS COSMOS 494057. The high-resolution ALMA [Cu] 
observations characterize the source as a prototypical high-z, tur- 
bulent, but rotating disk galaxy (Herrera-Camus et al. 2022; Par- 
lanti et al. 2023). These spatially resolved observations also re- 
veal the presence of a broader component in the spectrum of 
the [Cu] line emission which was interpreted as a neutral out- 
flow (Herrera-Camus et al. 2021). JWST observations HZA al- 
low us to investigate the warm ionized counterpart and compare 
its properties with those of the cold gas. 

This work is structured as follows: in Sect. 2 we describe the 
JWST and ALMA observations and data reduction, in Sect. 3 we 
describe the analysis of the integrated rest-frame optical spec- 
trum, in Sects. 4 and 6 we analyze the JWST and ALMA data at 
a spaxel level. In Sect. 5, we present the analysis of the spectral 
energy distribution (SED) and its results. Sect. 7 is dedicated to 
the analysis and the discussion of the multi-phase outflow prop- 
erties and their impact on the galaxy. In Sect. 8, we discuss other 
possible interpretations of the broad emission lines that we ob- 
serve, while in Sect. 9 we draw our conclusions. 

Throughout this work, we adopt the cosmological param- 
eters from Planck Collaboration et al. (2016): Ho = 67.7 
km s^! Mpc^!, Q,, = 0.307, and Q,= 0.691, giving 1”= 6.09 
kpc at z = 5.54. 


2. Observations 
2.1. JWST data 


HZA was observed with the NIRSpec instrument (Jakobsen et al. 
2022) onboard JWST as part of the GA-NIFS ? survey (program 
1217 - Integral Field Spectroscopy in COSMOS, PI: Nora Luet- 
zgendorf) on 23 April 2023. The galaxy was observed in the Inte- 
gral field spectroscopy (IFS) mode (Bóker et al. 2022) with both 
the G395H/F290LP and PRISM/CLEAR grating/filter combina- 
tion with an eight-dither position of the medium cycling pattern. 
The total integration time was ~ 5h for the high spectral resolu- 
tion observations (R ~ 1670—3600), which cover the wavelength 
range between 2.87—5.27 um, and ~ 1h for the prism observa- 
tions (R ~ 30 — 400), which span the wavelength range between 
0.6—5.3 um. 

We retrieved the raw data from the Mikulski Archive for 
Space Telescopes (MAST) archive, then we processed them with 
the JWST pipeline version 1.11.1 and Calibration Reference 
Data System (CRDS) context jwst_1097.pmap. We run the 
standard steps of the public JWST pipeline to produce the final 
cubes. In the pipeline process, we used some customized code 
to improve the quality of the final data (see Perna et al. 2023). 
In particular, during the first stage “calwebb_detector1” that ac- 
counts for detector level corrections, we also corrected for the 
1/f noise. At the end of the stage “calwebb_spec2” that cali- 
brates the data, we visually inspected the calibrated exposures, 
and we manually masked all regions affected by cosmic rays, 
persistences, and failed open shutters in the MSA mask. More- 
over, we removed the outliers by calculating the derivative of the 
count rate along the dispersion direction and removing all data 
where the measurement was above the 98th and 95th percentage 
of the distribution for R2700 and R100, respectively (D' Eugenio 


? https://ga-nifs.github.io 
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Fig. 1. Integrated [Om]45007À flux map of HZA from the JSWT NIR- 
Spec/IFS R2700 cube. We show the 5, 25, and 75c of [Om] as green 
contours. Red contours illustrate the [Cu] emission from the ALMA 
Briggs-weighted (robust = 0.5) cube, at 5, 15, and 250, respectively. 
HST/WFC3 F160W image is reported in black contours. The pink, 
white, and green ellipses represent the ALMA beam size, the HST 
F160W point spread function (PSF), and NIRSpec IFS PSF at 3.2 um 
estimated by D'Eugenio et al. (2023), respectively. 


et al. 2023). Then, we applied the third stage, “calwebb_spec3” 
to combine each exposure and create the final datacubes with a 
spaxel size of 0.05" by using the “drizzle” weighting. Finally, we 
subtracted the background emission from each cube by remov- 
ing the median spectrum computed in the target-free regions of 
the cubes. 


2.2. ALMA data 


The [Cu] raw data were retrieved from the ALMA archive 
(2012.1.00523.S PI: Capak; 2017.1.00428.L, PI: Le Fevre; 
2018.1.01605.S, PI: Herrera-Camus). The source was observed 
for ~ 20 min, ~ 30 min, and 4.7 h, respectively, in Band 7. We 
used CASA (McMullin et al. 2007) to combine these three ob- 
servations and calibrate the uv visibilities by using the script in- 
cluded in the dataset. The measurement sets were then processed 
with CRISTAL’s reduction pipeline (see Posses et al. 2024; Soli- 
mano et al. 2024; Villanueva et al. 2024; Herrera-Camus et al. 
in prep. for details). In brief, first, the continuum emission is 
subtracted in the visibility space using the task uvcontsub, the 
task tcleanis run with a Briggs weighting (robust parameter — 
0.5 and 2), hogbom deconvolver and a spaxel scale of 0.0445". 
The resulting datacubes have a beam size of 0.33" x 0.29", and 
0.45" x 0.41", respectively. The cube with Briggs robust param- 
eter — 2 (natural weighting) reaches a higher S/N ratio at the 
expense of the angular resolution, while the cube with Briggs 
robust parameter — 0.5 is a compromise between the natural and 
the uniform (high angular resolution but lower S/N) weighting, 
providing us a good S/N and a higher angular resolution with 
respect to the natural weighting cube. All cubes have a channel 
width of 40 km s7!. 


ALMA observations also allow us to detect the FIR dust con- 
tinuum emission at ~ 160 um in the galaxy rest frame. The con- 
tinuum image was created with the task tclean and a Briggs 
weighting with a robust parameter — 2 (natural weighting) to 
maximize the S/N. The final continuum image has a beam size 
of 0.46" x 0.41" and a sensitivity of 6 uJy/beam. 


2.3. Astrometric registration 


We verified that the astrometric coordinates of JWST are aligned 
with those from ALMA. Since many other NIRSpec IFS obser- 
vations showed an offset between the position indicated in the 
datacube header and the real position of the source on the sky 
(e.g., Arribas et al. 2023; Jones et al. 2023; Lamperti et al. 2024), 
we assessed whether the astrometry of the JWST cubes of HZ4 
was correct or not. 

First, we retrieved the F160W images of HST/WFC3 (PID: 
13641) from the MAST archive, with astrometry corrected a pos- 
teriori by using the GAIA eDR3 catalog with uncertainties of the 
order of 0.01". We then compared the image created by convolv- 
ing the R100 cube with the F160W filter response and the image 
from the F160W filter from HST. We fitted a 2D Gaussian profile 
to find the coordinates of the center both in the HST and JWST 
images. The mismatch between the 2 images corresponds to a 
0.10" shift. We then shifted the JWST R100 cube to match the 
F160W astrometry. We also ensured that the R2700 and R100 
cubes were consistently aligned by collapsing the [Om]45007À 
emission line visible in both cubes. After verifying the align- 
ment, we shifted the R2700 cube by the same amount as the 
R100 cube. We assumed that the ALMA astrometry is correct, 
considering that the ALMA positional accuracy varies from 5% 
to 20% of the beam depending on the S/N, which corresponds 
to an accuracy of 0.015"-— 0.060". The high S/N ensures that the 
ALMA astrometry is correct within 1 spaxel (i.e. 0.05"). Fig. 
1 shows the [Om] integrated emission along with the contours 
from the HST F160W image and the [Cr] integrated emission 
after the astrometric correction. 


3. Integrated analysis 


In this section, we describe the analysis of the R100 and R2700 
NIRSpec IFS spatially integrated spectra of HZ4. The spatially 
resolved analysis is reported in Sect. 4. 


3.1. Spectral fitting of the R2700 cube 


In the top panel of Fig. 2 we show the R2700 spectrum inte- 
grated over the regions of the cube with S/N > 3 in the wave- 
length range between 3.27 — 3.29 um, which encompasses the 
[Onr]45007À emission. The region from which we extracted the 
spectrum is highlighted in the inset panel as the red contours. 
The gray shaded area is the error on the cube computed as the 
quadratic sum of the error extension in the datacube rescaled to 
match the RMS in the line-free regions of the cube (the error was 
rescaled by a factor of ~ 3, see also Ubler et al. 2023). 

The G395H spectrum spans the wavelengths between 2.87 — 
5.14 um corresponding to 4388 — 7859 À rest-frame. In the inte- 
grated spectrum, we detect the brightest rest-frame optical emis- 
sion lines Ha, Hf, [Om]A45007,4959À, [Nr]446548,6584À, 
[$1]446717,6731À and He145875Å. 

We fitted the integrated spectrum by modeling it with the 
sum of a power-law continuum and two Gaussian components 
for each line, one narrow and one broad. The [S$1]446717, 6731 
and Her45875 emission lines were fitted only with one narrow 
line due to the low S/N. We found that one Gaussian profile 
is enough to reproduce the observed profile of these low S/N 
lines. We tied the kinematics of each set of line components 
(narrow and broad) to have the same velocity and line width. 
For the narrow component, we left the line width free to vary 
in the range 0 < c < 250 km s^!, while the broader compo- 
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Fig. 2. Upper panel: Integrated spectrum of the source extracted from the regions with S/N > 3 in the wavelength range 3.27 — 3.29 um 
encompassing the [Om]45007À emission. The extraction region is shown as the red image in the inset panel. Bottom panels: Zoom-in on the 
integrated spectrum (black) and our best fit (red). The narrow and broad components are shown as the green and blue dashed lines, respectively. 
The Hf and [Om] complex is shown in the left panel, and the Ha, [Nu], and [Su]446717, 6731 complex is shown in the right panel. The gray 


shaded areas mark the uncertainties on the data. 


nent was constrained to have a line width larger by at least 2096 
than the narrow one, and we imposed the amplitude of the broad 
component to be lower than 50% of the narrow one, following 
Carniani et al. (2023). The ratios of the intensities of the [Om] 
and [Nr] doublets have been furthermore constrained to atomic 
physics motivated ratios (see Osterbrock & Ferland 2006) of 
[Onr]45007/[0m]44959 = 2.98 and [Nu]46584/[Nn]46548 = 
2.94. 

We have also fitted the spectrum by adopting only one Gaus- 
sian component for each emission line, but this model has been 
revealed to be inadequate to reproduce the data. Indeed, the 
Bayesian information criteria (BIC) test? (Liddle 2007) highly 
favors the two components model with ABIC = BIC; — BIC, > 


3 By assuming a Gaussian noise, BIC = x? + kIn(N) where k is the 
number of free parameters of the fit, and N is the number of data points. 
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10 (Kass & Raftery 1995), where BIC; and BIC, are the val- 
ues of the BIC for the fit with one and two components for each 
emission line, respectively. 


In the bottom panels of Fig. 2 we show the best-fitting results 
and the residuals of the integrated spectrum analysis for the HB 
and [Om] complex on the left, and the Ha, [Nu], and [Su] com- 
plex on the right. The flux of each emission line and the FWHM 
of the narrow and broad components are reported in Table 1. 
This analysis allows us to compute the integrated properties of 
the source and compare it with all previous studies that could not 
resolve the galaxy morphology. From the narrow component, we 
estimate a redshift of 5.54455 + 0.00002 for the galaxy, which is 
in agreement with the redshift measured from the [Cn] emission 
(Béthermin et al. 2020; Herrera-Camus et al. 2021). On the other 
hand, a broader component (FHWM ~ 430 km s^!) is clearly 
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Table 1. R2700 line fluxes. 


Line Flux x10? [erg s lem ?] 
Hoa- Narrow 1617 +35 
Hf- Narrow 42] x 18 
[Om]45007À — Narrow 2268 + 42 
[Ni]46584À — Narrow 206 + 17 
Ha- Broad 986 + 75 
Hf- Broad 313 + 56 
[Om]45007À — Broad 1297 + 85 
[Ni]46584À — Broad 297 + 41 
He145875À 97 + 32 
[Su]46717A 122+9 
[S$1]46731À 107 x9 
Component FWHM [km s71] 
Narrow 180 x2 
Broad 430 + 10 


visible and needed to reproduce all the emission line profiles. 
This component is slightly blueshifted with respect to the narrow 
component by Vbroad — Vnarrow = —24 + 3 km s !. The FWHM of 
the broad component is comparable with the [Cn] one, observed 
by Herrera-Camus et al. (2021) in this galaxy, but we measure 
a blueshift, instead of the redshift measured by Herrera-Camus 
et al. (2021). We interpret this component as an outflow and we 
discuss its properties and possible interpretations in Sects. 7 and 
8. 

By assuming case B recombination, the theoretical ratio be- 
tween the Ha and Hf fluxes is Fyy/Fug = 2.86 (for a temper- 
ature T=10* K and a density n, = 10? cm^?, see Osterbrock & 
Ferland 2006, for details). Assuming the Calzetti et al. (2000) 
attenuation law (Ry = 4.05), we estimate the dust attenuation of 
the source from the measured Balmer decrement. We obtain a 
value of Ay = 1.0 + 0.2 for the narrow component nebular dust 
attenuation. By using the same method we also compute the ob- 
scuration of the outflow component, which is Ay — 037. 

To compute the star formation rate (SFR) of the galaxy, we 
use the calibration by Kennicutt & Evans (2012) that exploits the 
Ha emission line luminosity. We corrected the narrow Ha flux 
by the dust obscuration reported above, hence we compute an 
integrated SFR = TIE Mo yr !, which is consistent within the 
error with the one derived from the SED fitting by Faisst et al. 
(2020b) by exploiting optical and UV photometry, of 41733 Mo 
yr. 

By exploiting the ratio between the fluxes of [S1]46717À 
and [S1]46731À, we computed the electron density of the source 
Sanders et al. (2016) and we obtain a value of n, — 2707 cm. 
The estimated electron density is, within the large uncertainties, 
in agreement with the one obtained from [S11]446717, 6731 from 
a sample of SF galaxies at 2.7 « z « 6.3 (Reddy et al. 2023), and 
other high-z galaxies targeted in individual NIRSpec-IFS studies 
(Jones et al. 2024; Lamperti et al. 2024; Rodríguez Del Pino et al. 
2024). 


3.2. Spectral fitting of the R100 cube 


Following the same approach as for the high-resolution cube, we 
extract the spectrum from the prism cube from the same aper- 
ture as the R2700 one. We show the R100 integrated spectrum 
in Fig. 3. The prism spans the wavelength range between 0.6 
and 5.2 um allowing us to probe both the rest-frame UV and the 


Table 2. R100 line fluxes. 


Line Flux x10? [erg s lem ?] 
Ha 3032 + 110 

Hg 761 + 36 

Hy 325 + 35 
[Onr]45007À 3663 +71 
[Nu]46584À 200 + 91 
Her45875À 158 + 29 
[Nem]43870À 185 + 40 
[On]A443727,3729À 1153 + 90 
[S1]446717, 6731A 242 + 35 


optical region of the spectrum. In particular, we can now probe 
the rest-frame wavelengths « 4500 À, which are inaccessible in 
the R2700 cube. In the new range of wavelengths, we detect the 
[Or]443727,3729, [Nem]A3870 and Hy emission lines, as well 
as a blue UV continuum. We fit the integrated spectrum as a sum 
of a continuum emission and a set of emission lines with sin- 
gle Gaussian profiles, since the lower resolution of the spectrum 
does not allow us to distinguish between the outflow and the nar- 
row component, as we did for the R2700 cube. To account for 
the Balmer break, we modeled the continuum emission as two 
different power laws with a discontinuity at 3645À rest frame, 
then we convolved the spectrum with the line spread function of 
the prism (Jakobsen et al. 2022). The velocities of the lines were 
tied together, but the velocity dispersion was left free to vary in 
order to allow for the varying spectral resolution of the PRISM 
with wavelength. It is important to note that the lines are unre- 
solved in the prism spectrum, hence the line width only reflects 
the instrument's resolution. 

In addition, we set the flux ratio of the [Om] and the 
[Nn] doublets according to atomic physics prescriptions (see 
Sec 3.1). We use only one Gaussian component to fit the 
[$1]446717, 6731 and the [Or]443727,3729 doublets due to the 
low resolution, which does not allow us to deblend the contri- 
bution from each line of the doublet. Therefore, we obtain only 
the total flux of the doublet. Similarly, the [Nu] emission lines 
are blended with the Ha emission line, hence it is not possible to 
disentangle their contribution, and the Ha flux is not as reliable 
as the one estimated from the high-resolution spectrum. 

We report the line fluxes measured from the R100 spec- 
trum in Table 2, and note that within the errors, they agree with 
those reported for the R2700 cube (after summing the contri- 
bution from the narrow and broad components). This demon- 
strates that the flux calibration largely agrees between the two 
cubes. The only exception is the flux for the Ha and [Nri] emis- 
sion lines, which are blended together, however the summed 
flux agree within the uncertainties. From the expected line ra- 
tio of Fuy/Fug = 0.466, we estimated a dust extinction of 


Ay = Dus assuming that the [Om]44363À emission line flux, 
which is blended with the Hy in R100 observations, is zero. We 
find that this value agrees within the large uncertainties with the 
one computed from the Ha and Hf ratio from the grating. From 
now on we will use as fiducial value of Ay the one inferred from 
the R2700 grating. 

By exploiting the strong-line metallicity diagnostics based 
on rest-frame optical line ratios, we can calculate the gas 
phase metallicity of the galaxy. We adopt the calibrations 
from Curti et al. (2017, 2020), recently revised for high- 
z, low-metallicity galaxies (Curti et al. 2024; Laseter et al. 
2024). In particular, we exploited the diagnostics R3 = 
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Fig. 3. Prism spectrum integrated from the region marked by red contours in Fig. 2, and best-fit results. In black we show the data, while the gray 


shaded region represents the 10 error. In blue we show the best-fit emis 


[Om]A5007A/HB, R2 = [On]AA3727,3729/HB, R = 0.47 R2 + 
0.88 R3, O32 = [Omr]/[Or]443727,3729, N2 = [Nn]/Ha, and 
O3N2 = ((Om)]/H6)/([Nu]/Ha) to estimate the abundance of oxy- 
gen log(O/H). 

The [Or]443727,3729 emission line doublet is only covered 
by the R100 observations, which do not allow us to disentangle 
the outflow contribution. Hence, we estimated a global metal- 
licity by assuming by assuming the total line flux for each line 
(from the PRISM for [Or]443727,3729, narrow+broad from the 
grating for the other lines). We correct each line flux for the 
dust obscuration computed from the Balmer decrement com- 
puted as the ratio between the total fluxes (narrow+broad) of 
Ha and Hf in the grating (Ay = 0.7 + 0.3). We estimate a 
value of 12 + log(O/H) = 8.34 + 0.08 for the entire galaxy. 
By using only the R3, N2, and O3N2 calibrators (i.e. those 
provided by the grating spectrum) we are also able to estimate 
the galaxy versus outflow metallicity. We obtain an abundance 
of 12 + log(O/H)parow = 8.3 + 0.1 for the entire galaxy and 
12+log(O/H)outtow = 8.4 +0.1 for the outflow component. Both 
values are consistent with the value we infer by adding the cali- 
brators that exploit the [Om]AA3727,3729 lines. 

The value that we obtain is in agreement with metallic- 
ity obtained for other high-z galaxies with similar stellar mass 
dog(M,/Mo) = 97507 for HZ4, see Sect. 5) and merging 
systems (Nakajima et al. 2023; Curti et al. 2024; Sanders et al. 
2024; Venturi et al. 2024). Assuming a solar oxygen abundance 
of 12 + log(O/H) narrow = 8.69 (Asplund et al. 2021), we obtain a 
metallicity of Znarrow = 0.4 + 0.1 Zo and Zoutow = 0.5 + 0.1 Zo 
for the narrow and outflow components, respectively. 


4. Spatially resolved analysis 
4.1. Channel maps and continuum emission 


The high-resolution NIRSpec R2700 datacube shows that the 
galaxy is made up of different components at different veloci- 
ties. This can be seen from the ten channel maps of [Om] at dif- 
ferent velocities shown in Fig. 4. We set the zero velocity at the 
systemic redshift measured from the narrow component of the 
integrated R2700 spectrum. In each channel map, we show the 
contours at 5, 25, and 75o of the map at v = 0 km s^! to high- 
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sion lines, and in red the best-fit model. 


light the different morphologies of the [Om] emission at different 
velocities. 

The v = 0 map shows an elongated morphology, with a cen- 
tral peak that coincides with the peak shown in the UV observa- 
tions (see also Fig. 1) and an extended tail in the southern direc- 
tion. From the map at v = 60 km s^! it is evident that the elon- 
gation of the source is due to the presence of three components. 
In particular, a large southern clump is visible in the redshifted 
channels (v > 60 km s^!). We label this component as HZ4- 
South. The presence of these different components hints that the 
galaxy may be a galaxy merger or a clumpy galaxy, as observed 
for other high-z sources (Carniani et al. 2018; Arribas et al. 2023; 
Jones et al. 2023; Venturi et al. 2024). On the other hand, in the 
blueshifted channels, a component extending in the north direc- 
tion and peaking in the maps at v < —182km s^! slightly north- 
ern than the 750 contours peak at v = 0 becomes evident. This 
component hints at the outflow that we will investigate in the 
next Sections. 

We also looked at the continuum emission from these 
sources to identify these three different components. Fig. 5 
shows the continuum emission in three different wavelength 
ranges obtained from the NIRSpec R100 cube. In the left panel, 
we show the emission of the rest-frame UV continuum (1-2 
um observed wavelength or 1530-3058 À rest-frame). The UV 
emission morphology is similar to the one of the [Om] line 
(shown as the black contours); the location of the two peaks co- 
incides. Another source appears visible leftward of the target. 
We identify this source as a foreground source at z ~ 2.6 and 
we discuss in more detail its properties in Appendix A. From the 
rest-frame optical continuum in the observed range 3.3-4.2 um 
(5045-6422 A rest-frame) and 4.4—5.2 um (6730-7950 A rest- 
frame), shown in the central and right panels of Fig. 5, respec- 
tively. Here the southern, redshifted clump, already visible in the 
channel maps at v = 121 km s^, is now evident. This clump is 
very faint in [Om] emission (and any other emission line) but 
appears visible in continuum. The location of the central, small 
clump already identified in the channel maps at v = 60 km s^! 
is coincident with the location of the peak of the [Cu] emission 
(see Fig. 1) and the peak of the SFR discussed in Sect. 4. 

From the maps presented in this section, we can distinguish 
three different galaxy components. A northern component (HZA- 
North) which is very bright in rest-frame UV and optical con- 
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Fig. 4. Channel maps from the R2700 cube targeting the [Om]45007À emission line. Black lines are the contour of the 5, 25, and 75 times the 
RMS level computed in the line-free regions of the cube for the emission at v — 0. The velocity marks the central velocity of the channel. The 
black cross shows the location of the [Cu] emission peak (see Fig. 1). The white circles represent the position of the components of HZ4. 


tinuum and [Orr] emission, a central component (HZA-Central) 
which is bright in [Cir] emission, and a redshifted Southern com- 
ponent (HZ4-South), which is faint in [Om] emission but with an 
extended tail visible only in the rest-frame optical continuum. 


4.2. Spaxel-by-spaxel R2700 results 


The high S/N of the cube at the spaxel level allows us to perform 
a spatially resolved analysis of the galaxy. We fitted the cube at 
a spaxel-by-spaxel level with both a single Gaussian component 
for each emission line and two Gaussian components, following 
the same approach discussed in Sect. 3.1. We then selected for 
each spaxel the best fit by performing the BIC test and selecting 
the model with two components for each line only when BIC; — 
BIC, > 10, which represents positive evidence toward the model 
with two components ( e ). In the other cases, 
we selected the model with only one component, as the addition 
of the broader component was not significantly improving the 
quality of the fit. We excluded all the spaxels with S/N « 3 for 
each emission line. The S/N is computed as the ratio between the 
flux of the best-fit emission line and the resulting error from the 
fit. 

The results of the fitting (flux, velocity, and velocity disper- 
sion maps) for the narrow and outflow components are shown in 
Fig. 6. The galaxy morphology traced by the rest-frame optical 
emission lines differs from what is observed in the rest-frame 
UV continuum and the [Cri] FIR emission line, as already seen 
in the channel maps in Fig. 4. The narrow [Om] flux map shows 
a bright clump in the northern region (HZA-North). This clump 
is in the same location as the bright emission observed in the 
F160W filter from HST (see also Fig. 1), which is tracing rest- 
frame UV emission. Additionally, the [Our] map shows an ex- 
tended tail, in the southwest direction. The tail in [Om] shows 
the presence of a second faint component (HZ4-Central), which 
is instead more prominent and visible in the Ha map, and is also 
where the [Cu] peak. This suggests that the second component 
is likely tracing a region with higher dust obscuration, as we dis- 
cuss later in this section. 

The narrow component velocity field (top-central panel of 
Fig. 6) shows a different picture from what was seen in ALMA 


observations. Whereas the [Cu] kinematics is consistent with a 
rotating disk ( ), the ionized gas shows 
a more complex kinematic pattern, thanks also to the higher spa- 
tial resolution of the JWST data. In the northern region of the 
galaxy, we observe a relatively shallow velocity gradient rang- 
ing from 50 to -20 km s^!, while the southern region is red- 
shifted by ~ 100 — 130 km s^!. The clumpy morphology of the 
galaxy, combined with the observed abrupt velocity gradient in 
the southern part, suggests that the galaxy is undergoing a merger 
event rather than being a rotating disk as inferred from ALMA 
Observations ( ; : 
j. 

In the central row of Fig. 6 we show the outflow maps. In the 
outflow flux map, we also show as green lines the narrow com- 
ponent 5, 25, and 750 contours. The outflow morphology shows 
an elongated shape in the same direction as the narrow compo- 
nent emission. The brighter region is also in the northern part, 
but the peaks of the broad and the narrow emission are not in 
the same location. It is evident that the narrow component peaks 
~ 0.1” to the south of the outflow emission, which corresponds 
to a projected distance of 0.6 kpc. The velocity map shows that 
in the majority of the spaxels the outflow is blueshifted by up 
to -130 km s^! relative to the galaxy-integrated systemic red- 
shift computed from the narrow component. The velocity shift 
relative to the narrow component increases in the southern part 
of the outflow component. The southern part coincidentally has 
also the highest velocity dispersion. 

In Fig. B.1, we show the de-projected outflow velocity calcu- 
lated as Vout = |AVnarrow,broad| +20 out. The velocity reaches its peak 
of approximately 600 km s^! at the position of the HZ4-Central 
component. This location coincides with the highest Uspr, which 
is indicated by dashed contours in Fig. and will be discussed 
at the end of this section. The presence of the higher velocity 
at the location of the maximum SFR and the fact that the peak 
flux of the broad component does not coincide with the peak of 
the narrow one lead us to interpret this outflow as launched from 
the HZ4-Central component. Another possibility could be that 
we are observing multiple winds driven by the SF regions dis- 
tributed in the system. We will delve deeper into the geometry 
and the outflow properties in Sect. 
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Fig. 5. Continuum emission from the R100 JWST/NIRSpec datacube. From left to right the average continuum emission between 1—2 um (1530- 
3058 À rest-frame), 3.3-4.2 um (5045-6422 A rest-frame), 4.4-5.2 um (6730-7950 A rest-frame), respectively. The magenta circles represent the 
apertures from which we extract the spectrum for the North, Central, and Southern components for the analysis. The black contours represent the 
[Om] emission at 5, 25, and 75 c (see Fig. 1). The black cross shows the location of the [Cu] emission peak (see Fig. 1). The x and y axes represent 
the displacement in arcsec with respect to the galaxy center (RA = 09h58m28.5s, Dec= +02d03m06.7s, see Sect. 1). 
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Fig. 6. Results for the spaxel-by-spaxel fitting of the R2700 cube. In the top panels, we show from left to right the [Om]45007À flux, velocity and 
velocity dispersion maps of the narrow component, respectively. In the central panels, we show the same but for the [Om] outflow component. In 
the [Om] outflow flux map, we show the 5, 25, and 750 contours of the narrow [Om] flux in green. In the bottom panels, we show from left to 
right the flux maps of the Ho, Hf, and [Nr]46584À narrow components, respectively. The plus signs show the location of the three components 
identified in Figs. 4 and 5. 
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Fig. 7. Top panels: BPT diagram on the left and VO87 diagram on the right. As pink and turquoise crosses we show the location on the BPT 
diagram of the narrow and the outflow component for the integrated spectrum. As red and green points we report the location in the BPT diagram 
of each spaxel for the narrow and outflow component, respectively. The gray shaded region shows the location of SDSS galaxies, while blue 
diamonds show the location in this diagram of low-metallicity high-redshift (z ~ 5.5 — 9.5) star-forming galaxies from Cameron et al. (2023). As 
black dashed and dast-dotted lines, we report the demarcation line for the local star-forming galaxies and AGN samples by Kewley et al. (2001) 
and Kauffmann et al. (2003). As the dash-dotted black line we show the demarcation line between high-z low-metallicity galaxies and AGN by 
Scholtz et al. (2023). AGN are right-ward of the line, while on the left there is a mixed population of AGN and SF galaxies. Bottom panels: from 
left to right the logarithm of N2, R3 in each spaxel for the narrow component, and the logarithm of N2, R3 in each spaxel for the broad component. 


In Fig. 6 we also present the [Nu] flux map, which shows 
a different morphology from the flux maps of the other emis- 
sion lines. We can highlight the presence of two peaks in the 
[Nu] emission, one associated with the northern component, and 
one associated with the central one. Differently from the Ha and 
[Om] the brightest [Nm] emission arises from the central region. 
These differences can be a sign of a different metallicity between 
the two galaxies or a different ionization mechanism, with the 
presence of a harder radiation field arising in the center. 


To study this latter possibility we exploit the Baldwin- 
Phillips-Terlevich (BPT; Baldwin et al. 1981) and the VO87 
(Veilleux & Osterbrock 1987) diagnostic diagrams, which are 
sensitive to the ionizing source in a galaxy, between star forma- 
tion in HII regions and AGN photons. In particular, the first one 
is based on the line ratios of [Nu]/Ha and [Om]/Hf while the 
latter uses the [S$1]446717, 6731/Ha and [Om]/H£. 


In the bottom-left panel of Fig. 7 we report the spatially in- 
tegrated measurement of the galaxy in the BPT diagram for the 
narrow and the outflow component, together with the location of 
each spaxel. We observe that the integrated emission and the ma- 
jority of the spaxels, of the narrow component, lie on the Kewley 


et al. (2001) demarcation line, together with the low-metallicity 
z —5.5-9.5 galaxies from Cameron et al. (2023). On the other 
hand, the majority of the outflow component points lie in the 
composite region of the BPT between the demarcation lines of 
Kewley et al. (2001) and Kauffmann et al. (2003) that are based 
on the SDSS galaxy sample (gray shaded region). Low-redshift 
observations have shown that in the presence of shocks induced 
by starburst-driven winds, or by tidal flows due to a merger, the 
line ratios can shift toward the composite region even in the ab- 
sence of an AGN (Monreal-Ibero et al. 2010; Rich et al. 2011, 
2015; Medling et al. 2018). The demarcation lines used in the 
BPT diagram for low-z sources have proved to be unreliable at 
high-z (z > 3). Due to the increase of the ionization parameter 
and the decrease of the metallicity at high redshift, both galaxies 
and AGN may end up in the same location in the BPT diagram 
(Nakajima & Maiolino 2022; Übler et al. 2023; Harikane et al. 
2023; Maiolino et al. 2023; Scholtz et al. 2023). We then used 
the new demarcation line by Scholtz et al. (2023), shown as a 
dash-dotted line in the plot. According to this new line, there 
are no points that lie in the AGN region, but all points lie in the 
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Fig. 8. Spatially resolved properties of the narrow component. From left to right: gas-phase metallicity (12 + log(O/H)), dust extinction (Ay), and 
dust-corrected star formation rate surface density (Zsrg) maps. The plus signs show the location of the three components identified in Figs. 4 and 
5. The gray contours represent the 5, 25, and 750 of the narrow [Om] flux from Fig. 6. 


region in which SF galaxies and AGN may have the same line 
ratios. 

In the bottom-right panel, we show where the narrow compo- 
nent lies in the VO87 diagram, together with the sample of high-z 
galaxies from Cameron et al. (2023) and the local galaxies from 
SDSS. The point lies on the demarcation lines by Kewley et al. 
(2001) and Scholtz et al. (2023), therefore we cannot constrain 
the ionization source. 

Unfortunately, for the wavelength range and spectral resolu- 
tion available, there are no other diagnostic diagrams that have 
been proven to be reliable at high redshift (Scholtz et al. 2023; 
Mazzolari et al. 2024). The absence of Hen and other high- 
ionization lines rule out the presence of a strong, energetic AGN 
and favor a SF interpretation, although we cannot completely 
rule out a contribution to ionization from a weak, subdominant 
AGN. 

We display the maps of the log([Nn]/Ho) and log([Omr]/H) 
ratios for the narrow and outflow components in the top row of 
Fig. 7. For the narrow component, a gradient in the values of 
log([Nn]/Ho) is evident, with the ratio increasing from north to 
south, while the values of log([Onr]/Hf) are almost constant. The 
outflow component, on the other hand, shows a gradient also 
in the log({Om]/Hf) values, varying from 0.2 to 0.7, while the 
log([Nn]/Ho) values range from -0.8 to —0.2. Since the AGN 
scenario is disfavored, the observed gradients of the line ratio 
can be interpreted as a metallicity gradient between the southern 
and the northern clump with the northern clump being poorer in 
metals compared to the southern one, given that the [Nu]/Ha ra- 
tio decreases as the metallicity decreases (Pettini & Pagel 2004; 
Curti et al. 2017, 2020). Metallicity gradients have been ob- 
served at high redshift, especially in clumpy or merging sys- 
tems (Arribas et al. 2023; Rodríguez Del Pino et al. 2024; Ven- 
turi et al. 2024). The study of the metallicity gradients for this 
galaxy will be presented in Curti et al. (in prep.), but here we 
show the metallicity maps for the narrow component in the left 
panel of Fig. 8. We exploit the N2, R3, and O3N2 line ratios 
and we adopt the calibrations from Curti et al. (2017) to derive 
the O/H abundances. The estimated abundances vary between 
8.2 « 12 + log(O/H) « 8.5. The smallest metallicity values are 
near the center of the northern and central components, while for 
the southern component, the lack of detection of [Nu] emission 
does not allow us to constrain the metallicity. 

Having excluded the presence of an AGN, we can also es- 
timate the SFR surface density (Zsrg) at a spaxel level by ex- 
ploiting the Ha dust-corrected luminosity. First, we estimated 
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the dust attenuation at a spaxel level thanks to our detection 
of Ha and Hf emission lines. By assuming case B recombi- 
nation, the theoretical ratio between the Ha and Hf flux is 
Fgo/Fgg = 2.86 (Osterbrock & Ferland 2006). We measured a 
spaxel-by-spaxel ratio ranging from 2.86 to 5.60. Assuming the 
Calzetti et al. (2000) attenuation law, we estimate a dust attenu- 
ation of 0 € Ay x 2.5. In the central panel of Fig. 8 we show the 
map of dust attenuation across the galaxy. 

Interestingly, the northern region, where the line emission 
and the UV continuum peak, shows little to no dust extinction, 
similar to what is observed in the southern component. How- 
ever, there is a zone of high extinction between these two regions 
where the dust attenuation is reaching values of Ay ~ 2.5 which 
is cospatial with the central component. 

By correcting the Ha flux for the dust extinction and using 
the calibration from Kennicutt & Evans (2012) to get the SFR 
from Ha luminosity, we computed the resolved SFR surface den- 
sity map, which is shown in the right panel of Fig. 8. This shows 
that there is a region of intense and obscured ongoing star for- 
mation in the galaxy, which is coincident with the location of the 
central component. By comparing with the [Cu] flux contours 
shown in Fig. 1, it can be seen that this region is coincident with 
the peak of the [Cri] emission, which indeed is a tracer of star 
formation unaffected by the presence of dust (De Looze et al. 
2011; Pineda et al. 2014; Herrera-Camus et al. 2015). 


5. SED fitting 


The prism data allow us to perform a spectral energy distribu- 
tion (SED) fitting to study the galaxy's stellar population and 
infer key galaxy properties such as the stellar mass and the star 
formation history. We used the SED fitting code Bagpipes (Car- 
nall et al. 2018). We adopted the stellar population SED model 
from Bruzual & Charlot (2003), and we used Cloudy (Chatzikos 
et al. 2023) to obtain the models for the nebular emission with a 
range in the ionization parameter priors of —3 < log U < 0. We 
assumed a Calzetti et al. (2000) dust attenuation law, with Ay 
free to vary between 0 and 4. We assumed the same extinction 
for the nebular and stellar components. We left the metallicity 
free to vary within the range 0.1 — 1 Zo. For the star formation 
history (SFH), we adopt a so-called non-parametric SFH, with 
continuity priors (Leja et al. 2019) with 10 bins logarithmically 
spaced from 1 Gyr (z = 5.5, the redshift of the source) up to 200 
Myr after the Big Bang. 
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Fig. 9. On the left: R100 spectra in black and SED best-fitting result in red extracted from the entire system, the northern clump, the central clump, 
and the southern clump, from top to bottom, respectively. On the right: best-fit star formation history. 


First, we performed the SED fitting on the integrated R100 
spectrum of the regions in which the S/N on the [Om] emis- 
sion line is larger than 3 (see red region in the inset panel in 
Fig. 2). We show the spectrum and the best-fit model in the 
top-left panel of Fig. 9, while we show the SFH in the top- 
right panel. The properties resulting from the SED fitting are 
reported in Table 3. We infer a total stellar mass of the system 
of log(M,/Mo) = 975 which is smaller than the previous 
estimate from the literature of log(M,/Mo) = 10.15 + 0.15 of 
Faisst et al. (2020b), but in agreement with the one from Capak 
et al. (2015) of log( M, /Mo) = 9.67 + 0.21. 


The SFH for the whole HZA system shows an almost con- 
stant SFR of a few Mo yr! up to the last 20 Myr, where we see 
a rapid increase reaching ~ 100 Mo yr^! that could be triggered 
by the merger (Pearson et al. 2019). Indeed, the integrated spec- 
trum shows the Balmer break at ~ 2 um, which is indicative of 
an older stellar population, while the luminous Ha is indicative 
of SFR in the last 10 Myr. 


If we extract the spectra for the northern, central, and south- 
ern components (see apertures in Fig. 5) we see that the older 
population is dominant in the southern component, with a well- 
defined Balmer break, faint continuum UV emission, and faint 
emission lines. On the other hand, the spectra extracted in the 


central and northern regions show a bright UV continuum. This 
difference in the spectra indicates a different stellar population 
that reflects different formation histories of the northern, central, 
and southern components (see Fig. 9, right panels). The south- 
ern component started to form stars only ~ 200 Myr after the Big 
Bang and has continued to grow with a slowly increasing SFR, 
reaching values of SFR = 8 + 2 Mo yr’! at the time of the ob- 
servation. From the SED fitting we obtain a mass-weighted age 
of the southern galaxy of 2s Myr. On the other hand, the 
northern galaxy shows a much quicker rise in the SFR, having 
assembled most of its mass in the last 30 Myr. The central com- 
ponent has a rising SFH very similar to the one of the northern 
one. Both the northern and the central components have a very 
young stellar population, with a mass-weighted age of « 100 
Myr. 


The differences in the SFH and the stellar populations sup- 
port our interpretation that this galaxy is formed by merging sys- 
tems, rather than three star-forming clumps. Cosmological sim- 
ulations from Nakazato et al. (2024) show that a merger event 
can induce the formation of clumps which lead to a bursty star 
formation with timescales of 10-30 Myr, consistent with what 
we observe in the northern and central clumps. 
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Fig. 10. Example of best-fit results for the fitting of the [Cr] emission line in one spaxel. On the left, we show the best-fit results by using two 
components, on the right the best-fit results by using one component. In the bottom panels, we show the residuals of the fitting procedure. The 


dotted lines show the RMS computed in the line-free regions of the cube. 


Table 3. Galaxy properties inferred from the SED fitting. 


log( M, /Mc) Ay T/Myr 
(1) Q) (3) 
HZ4 System | 9.75077 [| 0.39 + 0.04 | 357782 
HZ4 North | 8.90x0.05 | 0.28+0.03 | 27:8 
HZ4 Center | 8.69+0.06 | 0.25 + 0.04 50710 
HZ4 South | 9.24+0.08 | 0.71 +0.13 322788. 


Notes. (1) Stellar mass; (2) Dust attenuation; (3) Mass-weighted stellar 
age. 


6. ALMA analysis 


In this section, we describe the analysis of the ALMA obser- 
vations of the [Cir] emission line and the dust continuum. The 
ALMA [Cir] cube was analyzed with the same method as the 
NIRSpec R2700 cube, to detect and possibly spatially resolve 
the neutral outflow. Since the properties of the outflows inte- 
grated in one beam are already discussed in Herrera-Camus et al. 
(2021), here we focus only on the spaxel-by-spaxel analysis. To 
study the outflow component we exploit the datacube which was 
reduced with natural weighting, which has a higher S/N ratio 
at the expense of angular resolution with respect to the Briggs 
datacube. We fitted each spaxel having a S/N of the [Cm] emis- 
sion greater than 3 with both a single and two Gaussian lines 
and a constant continuum. The continuum emission was already 
removed from the cube during the data reduction, but we added 
this component in the fitting to furthermore remove any residual 
continuum emission that was not accounted for. Using a similar 
method as described for the JWST datacube, for each spaxel we 
selected the 2-component model only when BIC; — BIC, > 5 
and the S/N of both the narrow and the broad component were 
greater than 3, otherwise the 1-component model was adopted. 
We note that in this case, we use a slightly lower threshold for 
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the BIC selection due to the lower S/N of the line, but ABIC » 5 
still gives positive evidence in favor of the two Gaussian model. 
We applied the same procedure to the Briggs (robust — 0.5) dat- 
acube. Unfortunately, the lower S/N of this latter datacube does 
not allow us to select the outflow component in any spaxel due 
to its S/N falling below our threshold of 3, even if the ABIC 
was larger than 5 in some spaxels. But a broad component is 
still detected integrating into circular apertures as also shown by 
Herrera-Camus et al. (2021). 


6.1. [Cu] spaxel-by-spaxel results 


Here we present the results of the spaxel-by-spaxel analysis of 
the ALMA data targeting the [Cir] emission line. In Fig. 10 we 
show the best-fit results for the [Cm] emission line in one spaxel. 
We present the best fit results by using one component to fit the 
emission line in the right panel, and by using 2 components in 
the left panel. The 1-component fit shows an excess around a 
velocity of ~ 200 km/s. These residuals completely vanish by 
adding the second component. The 2-component model is fa- 
vored by the BIC test, with a ABIC — 14. We perform this test 
on each spaxel and we end up selecting the outflow component 
for 115 spaxels. 

In Fig. 11 we show the results for the spaxel-by-spaxel fitting 
of the [Cm] data. In the first two rows we show moment maps of 
the narrow and outflow components from the natural weighting 
reduced cube, while in the panels in the bottom row, we show the 
moment maps from the datacube reduced with Briggs weighting, 
robust-0.5. 

The cube reduced with natural weighting has a larger beam, 
but a higher S/N on a spaxel level with respect to the Briggs, ro- 
bust=0.5 one. The higher S/N allows us to select an outflow com- 
ponent, that we show in the central row. This component is red- 
shifted with respect to the host galaxy. This redshifted outflow 
is what was already observed in circular apertures by Herrera- 
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Fig. 11. [Cu] moment maps derived from a spaxel-by-spaxel Gaussian fitting. From left to right we present the maps of the amplitude of the 
Gaussian, the velocity, and the velocity dispersion, respectively. From top to bottom, we show the narrow and the outflow components from the 
natural cube and the total emission from the 1-component fit from the Briggs robust=0.5 cube. We note that no outflow components were selected 
from the Briggs cube because they fall below the 3c threshold. In the left panels we show the ALMA beam as a green ellipse and the narrow and 


outflow [Orr] flux from Fig. 6 as lime and red contours, respectively. 


Camus et al. (2021), similarly elongated along the northwest di- 
rection. In contrast to the [Om] outflow, it is redshifted with re- 
spect to the host galaxy (Av ~ 35 km s^!), but the projection 
is co-spatial with the [Om] outflow (shown as the red contours). 
The broad [Cu] component appears extended toward the filamen- 
tary structure (feature in the North-West at RA ~ 9h58m28.47s, 
Dec ~ 2d0.3m07.0s) that we observe in the narrow component 
natural flux map and in the map of the total profile. The moment 
maps of the narrow component show an elongated [Cri] emis- 
sion due to the presence of two separate spatial components, one 
peaking where the [Om] peaks (the north component), one at the 
location of the central component, where the integrated [Cu] map 
from Fig. 1 peaks. The smeared velocity field resembles that of 
a rotating galaxy as already observed in other studies (Herrera- 
Camus et al. 2022; Parlanti et al. 2023). 


In the bottom panels of Fig. 11 we show the maps of the nar- 
row component obtained from the fitting of the Briggs cube. The 
Briggs robust = 0.5 cube unambiguously confirms the presence 
of two distinct clumps (northern and central component) in the 


amplitude map, through the higher spatial resolution. The peak 
[Cu] flux is in the location of the central obscured component. 
This is in agreement with the [Cn] being a tracer of star forma- 
tion as the maximum luminosity arises from the regions with 
higher star formation (see also Sect. 5). 


6.2. Dust continuum emission 


The ALMA observations also allow us to study the resolved dust 
continuum emission at ~ 160 um in the galaxy rest frame. While 
the optical and UV rest-frame emission allows us to probe di- 
rectly the light emitted from the young stars, the FIR continuum 
emission traces the thermal radiation emitted from the dust that 
is heated by the stars. As the UV emission is absorbed by the 
dust, we can probe the obscured SF in the system by analyzing 
the dust emission. 

In Fig. 12, we report the dust emission contours at 3, 6 and 8 
c on top of the Ay map. The dust emission morphology follows 
the elongation of the source, already seen in [Om] and [Cr] emis- 
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sion. The dust continuum peaks between the northern and the 
central components, where the Ay measured through the Balmer 
decrement is 0 instead. The beam size (0.46" x 0.42") is larger 
than the separation between the two components. The difference 
between the Ay map and the dust continuum emission can be 
ascribed to a difference in the column density and the dust tem- 
perature across the field of view. The total flux of the source at 
160 um is S160, = 134 + 33 uJy, comparable with the values 
reported in Herrera-Camus et al. (2021) and Mitsuhashi et al. 
(2023) with a similar resolution and Faisst et al. (2020a) with 
lower resolution. 


The observations of the FIR continuum, together with the 
UV continuum and Ha emission line, allow us to compare the 
SFR computed with different methods. Before the advent of 
JWST, the total SFR of high-z galaxies, due to the lack of op- 
tical rest-frame emission observations, was estimated by sum- 
ming the unobscured SFR estimated from UV continuum emis- 
sion, and the obscured one estimated from the FIR continuum. 
The three, low-resolution (~1’’), dust continuum detections in 
ALMA bands 6, 7, and 8 of HZ4 allow to perform a FIR SED 
fitting to estimate the obscured SFR. Faisst et al. (2020a) report 
SFRuy = 29 +4 Mo yr“! exploiting the UV HST observations, 
and through SED fitting they estimate a SFRir = TI Mo 
si 

Despite three FIR detections, the SFRyr has uncertainties of 
^ ] dex, principally driven by the uncertainties in determining 
the dust temperature (Tp = ad K; Faisst et al. 20202) due 
to the lack of observations near the IR peak, which can strongly 
constrain the temperature. 


Thanks to the JWST observations, we were able to determine 
the galaxy SFR thanks to the observed Ha luminosity, corrected 
for the dust obscuration according to the Balmer decrement. As- 
suming that the dust-corrected Ha luminosity is able to recover 
the entire intrinsic SFR ongoing in the system, we obtain that 
the total SFR of the system is SFRu, = 77*12 Mo yr. Of this, 
29 +4 Mo yr^! is unobscured SF, and the remaining is obscured 
SF of SFRobscurea = 487 Mo ve. By using the relation be- 
tween the obscured SFR and Lir by Kennicutt & Evans (2012) 
log (SFR/[Mo yr7!]) = log(Lyjr/[erg s^! ]) — 43.41, we infer an 
IR luminosity of log(Lig/Lco) = 11.5 + 0.2. The IR luminosity 
obtained is comparable within the errors to the value obtained in 
other works from the SED fitting (Faisst et al. 2020a; Fudamoto 
et al. 2023). 

Then, we can now perform a SED fitting for the three FIR 
detections (we adopt the values from Faisst et al. 2020a for the 
band 6 and 8 fluxes and the value reported above for the band 
7 detection), but constraining the Lip from the obscured SFR 
inferred above. We assumed that the IR and Ha trace SF over 
a similar timescale, a dust emissivity index of B = 2, and we 
performed a fit assuming a greybody curve with a single temper- 
ature. With these assumptions, we obtain a dust temperature of 
Tp = 70 + 12 K and a dust mass of log(Mp/Mo) = 6.62 + 0.17. 
The dust temperature is in agreement with the values of Faisst 
et al. (2020a) and Fudamoto et al. (2023), but with smaller un- 
certainties. 


7. Outflow properties 


In this section, we present the properties of the ionized and neu- 
tral atomic outflows detected in JWST and ALMA data. We in- 
vestigate the origin of the outflows and their launching mecha- 
nisms. Moreover, we present and discuss the outflow history. We 
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Fig. 12. Ay map from narrow component (same as in Fig. 8) and dust 
continuum emission contours at 160 um rest-frame. In black we report 
the dust continuum 3, 6, and 8c contours. The black plus signs represent 
the location of the three components identified from the [Om] channel 
maps and in optical continuum emission. The purple ellipse represents 
the beam size of the dust continuum emission. 


consider the possibility that these broad components may trace 
other non-circular motions in Sect. 8. 


7.1. lonized outflow properties 


We first computed the integrated outflow properties for the ion- 
ized gas phase. By using the results from the modeling of the 
integrated spectrum reported in Table 1, we can compute the 
ionized gas mass from the flux of the Ha broad emission line 
component following Cresci et al. (2023), as 


(1) 


L a. 1 -3 
Mionou = 3.2 X 10° | Ha, outflow Jmm, Mo. 


10% erg s7! Nne 


and the broad [Om] emission following Carniani et al. (2023), as 


M; -0.8 x 108 | Liorir]outfiow | (= x ( Zs 
1on,out — A 


Mo, 2 
10% erg s-! Z ) S Q) 


Ne 


where Lye, outfow and Lrom,outflow are the luminosity of the 
Ha and [Om] broad components corrected for the dust extinc- 
tion, n, the outflow electron density and Z the metallicity of 
the outflowing gas. We assume n, = 200 cm^?, which is con- 
sistent with the value derived from the [Sr]446717, 6731 ratio 
for the ISM of the galaxy with the integrated fit, and we adopted 
an outflow metallicity of Z = 0.5 Zo (see Sect. 3.2). By using 
the value of Ha and [Om] fluxes reported in Table 1, we cor- 


rected for the outflow extinction of Ay = 037. and we estimate 


an ionized outflow mass of Mionou = (1.82 5) x 107 M, and 
Mionout = (25:3) x 107 Me for Ha and [Orr], respectively. 

For the outflow velocity voy, we adopted the following rela- 
tion from Rupke & Veilleux (2013) and Fiore et al. (2017): 
(3) 


-1 
Vout — [AVnarrow,broadl *20,4-389r12kms , 
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where AVnarrow broad 18 the velocity shift between the peak of the 
broad and the narrow components and Cout is the velocity disper- 
sion of the broad component deconvolved for the instrumental 
broadening. With this assumption, we are accounting for the un- 
known geometry of the outflow, under the assumption that only 
the tail of the broad wings trace the outflowing gas directed in (or 
moving along) our line of sight, and thus its intrinsic velocity. 
The mass outflow rate is computed by assuming that the mass 
outflow rate is constant with time (e.g., Lutz et al. 2020) 


; M, 
Masc (4) 
Rout 
where Rout is the extension of the galactic outflow. Based on the 
spaxel-by-spaxel fit we have inferred a flux-weighted radius of 
Rou ~ 2 kpe (see Fig. 6) that yields Mourne = 205 Mo yr! 


6 
and Mow, (om = cu Mo yr t. 


7.2. Neutral outflow properties 


With the new analysis of the ALMA data, we can also re- 
estimate the neutral atomic outflow properties. We computed the 
neutral atomic gas mass in the outflow by using the following 
equation: 


My = kicin(T, ng, Z) X Licin (5) 
where kicn is a conversion factor that depends on the temper- 
ature (T), the neutral gas density (ny), and the gas metallicity 
(Z), and Lycrn is the luminosity of the observed [Cu] outflow. 
Since the outflow temperature and density are unknown we con- 
servatively compute the value of ktc assuming T = 10^ K and 
ng = 10* cm^, which correspond to the maximal excitation 
case. With this assumption, the mass obtained represents a lower 
limit on the neutral gas mass traced by the [Cn] emission line. 

Assuming a oxygen abundance of 12+log(O/H) = 8.34 (Z = 
0.45 Zo, see Sect. 3.2), the case for maximal excitation gives 
kicm = 4.88 Mo/Lo (Herrera-Camus et al. 2021). The mini- 
mum neutral mass in the outflow is then My out = 2.69 x 105 Mo. 
We compute the flux-weighted radius of the outflow from the 
map shown in Fig. 11, which gives Ruou ~ 1.3 kpc, and the 
outflow velocity computed as in Eq. 3 is vgou ~ 260 km s 
Finally, we compute a lower limit on the mass outflow rate of 
Myou ~ 121 Me yr^!. The value that we obtain is higher than 
the lower limit reported in Herrera-Camus et al. (2021) of ~ 34 
Mo yr`!, but this is principally due to the difference in the kjc;jj 
(Herrera-Camus et al. 2021 assumed Kc; = 1.5 Mo/Lo, which 
correspond to maximal excitation assuming with solar metallic- 
ity). 

In HZA, the neutral atomic mass outflow rate 1s more than 
one order of magnitude higher than the ionized mass outflow 
rate. A similar behavior is observed both at low and high red- 
shifts and for AGN and SF galaxies (Herrera-Camus et al. 2020; 
Fluetsch et al. 2021; Avery et al. 2022; Cresci et al. 2023; 
D' Eugenio et al. 2023; Belli et al. 2023). Analyzing a sample 
of local ULIRGs, Fluetsch et al. (2021) shows that the amount 
of mass expelled by the galaxy in the ionized phase is usually 
negligible with respect to the neutral and the molecular phases, 
the latter of which accounts for the majority of the outflow mass. 
We compare our results with others from the literature in Fig. 13, 
where we show that both AGN and SF galaxies at low and high 
redshift lie below the 1:1 relation between the ionized mass out- 
flow rate and the neutral one. 
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Fig. 13. Comparison between the ionized and neutral mass outflow rate 
color-coded by redshift. The pink star marks the measurements for HZ4 
from this work. The turquoise circles are the local ULIRGs studied 
in Fluetsch et al. (2021). The triangle, pentagon, diamond, and square 
display the results obtained for local and high-z AGN by Perna et al. 
(2019); Cresci et al. (2023); Belli et al. (2023); D'Eugenio et al. (2023), 
respectively. The black dashed line represents the 1:1 relation. 


7.3. Outflow strength and effect on the host galaxy 


To understand whether these outflows are able to leave the 
galaxy, and hence remove gas from the system and enrich the 
CGM and the IGM, we estimate the escape velocity. We fol- 
lowed the same approach as in Carniani et al. (2023). We mod- 
eled the galaxy gravitational potential as the sum of an exponen- 
tial disk profile potential for the baryonic component with mass 
M, + My (the dust mass, of ~ 108 Mo, is negligible; Pozzi 
et al. 2021) and a Navarro-Frenk-White (NFW; Navarro et al. 
1996) profile for the potential of the dark matter halo. 

We considered as M, the stellar mass of the whole system 
reported in Table 3. For the gas mass, we convert the [Cr] lu- 
minosity as Mgas = @jcn X Lie. We adopt atem = 31 Mo/Lo 
(Zanella et al. 2018). To compute the mass of the dark matter 
halo, we adopted the stellar-to-halo mass relation from Moster 
et al. (2013), and we compute the concentration parameter for 
the NFW profile from the relation between the mass of the dark 
matter halo and the concentration estimated from Dutton & Mac- 
ció (2014). 

We find that the velocity necessary for the gas to leave the 
galaxy and reach infinity is Vese ~ 800 km s^!. This velocity 
is higher than the estimated outflow velocity of the ionized and 
neutral gas, hence these outflows are leaving the launching site, 
but eventually, the gas may fall back on the galaxy. We note 
that these calculations assume that the outflow moves ballisti- 
cally (i.e., any other slowing mechanisms, such as the drag due 
to the gas encountered in the way, are not taken into account) 
and it is launched from the center of these idealized profiles. 

To understand the impact that the outflow has on the re- 
moval of gas from the galaxy we estimate the mass loading factor 
n= Mout/SFR. For the integrated values we obtain Nionizea ~ 3%, 
while reuraj ~ 150%. The ionized value is smaller than what 
found for other high-z, but less massive galaxies (7ionizea ~ 20096 
for galaxies with log(M,./Mo) ~ 7.5 — 8.5; Carniani et al. 2023). 
Both cosmological simulations and observations predict that the 
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mass loading factor should decrease with increasing stellar mass 
for star-formation driven outflows, with mass loading factors 
reaching less than unity at a mass of around log(M, /Mo) ~ 10, 
consistent with what we observe (Muratov et al. 2015; Nelson 
et al. 2019; Pandya et al. 2021; Carniani et al. 2023). 


7.4. Outflow driver and launching region 


To investigate the outflow driver we first computed the energetics 
of the outflow. We compute the kinetic and momentum rates as 


; 1 
Eo = 5 


; (6 


y 2 
MoutVout 


(7) 


For the ionized outflow, we DUAE Eoi He = 1.16 x 10%! 
ergs, pm jour = 1.57 x 10^! erg s7 ', and Poat He = 5.9 x 103 
gcms pP jon] = 8.1 x 10"? g cm s™ 2 with uncertainties of 
the order of 0.5 dex. 

Even though we cannot completely rule out the possible 
presence of an AGN (see Sect. 3.1), the outflow velocity and the 
outflow energetics are small in comparison to what we expect 
from AGN-driven winds (Fiore et al. 2017). On the other hand, 
the outflow can be driven by the stellar feedback in a particular 
combination of winds from massive stars, and type 2 supernova 
(SN) explosions that release energy and momentum into the ISM 
(Veilleux et al. 2005) as explained below. 

The SFR threshold for driving an outflow is considered to be 
around Xsgg ~ 1 Mo yr! kpc^? (Newman et al. 2012; Davies 
et al. 2019). The star-forming obscured clump (HZ4-C) has a 
SFR of ~ 20 Mo yr! in a region of radius ~ 0.9 kpc, resulting 
in a Xsgg ~ 4 Mo yr! kpc? which is high enough to launch the 
ionized outflow. 

If the outflow is driven by the stellar feedback from the burst 
in the central dusty region, we can compare the energy rate of 
the SN exploding in that region with the kinetic rate of the out- 
flow. We can compute the energy rate driven by SN explosions 
as Esnn = €EsniR, where e is the supernova efficiency, that is 
the efficiency of transferring kinetic energy to the ISM, Eswr is 
the total energy released in a SN explosion (10°! erg), and R 
is the supernova rate. Assuming a SFR of ~ 20 Mo yr !, for a 
Kroupa IMF, the rate of SN explosions is related to the SFR as 
R ~ 0.01 yr! x SFR / (Mo yr^!) (Mo et al. 2010). If we assume 
a SN efficiency of € ~ 0.1, we obtain Esnn = 7.3 x 10^! erg s " 
which is higher than the outflow kinetic rate (Eou = 1.6 x 107! 
erg s x Following Heckman et al. (2015), the momentum rate 
injected by a starburst in the ISM is Psp = 10% TSFR / (Mo yr ) 
g cm s^?. In the case of the central clump, we obtain Psg ~ 1035 
g cm s?. Both the momentum and energy released by star for- 
mation are higher than those inferred for the outflow, thus the 
SF feedback in the central, obscured clump can drive the ionized 
outflow that we observe in the rest-frame optical emission lines. 

Moreover, we observe that the outflow has higher Av and 
FWHM near the region of higher SFR (see Fig. 6 and 8), which 
corresponds to higher outflow velocity corrected for the un- 
known geometry. This is what we expect if that is the launching 
region and the outflow moves ballistically and slows down due 
to gravity, and possibly due to drag by gas encountered. 

Since the position of the peak of the outflow flux and the 
narrow [Om] flux do not match, we can infer that the outflow 
peak is not an artifact due to the higher S/N ratio of the emission 
line, but rather a real higher flux of the outflow in that region. 


Pout = MoutVout 
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Fig. 14. Ionized mass outflow rate map from the dust-corrected broad 
Ha emission for each radial shell. 


The interpretation of these data is that the dust-obscured clump 
launches the ionized outflow which moves toward the north in 
the foreground of the galaxy with respect to our line of sight. A 
schematic representation of this interpretation is presented in the 
top-left panel of Fig. 17. 


7.5. lonized outflow history 


Speculating that the ionized outflow is being launched by the 
central, obscured, highly star-forming component, we extracted 
integrated spectra from 6 concentric annuli centered on that 
clump and spaced in radius by 0.1", which corresponds to 609 
pe at z = 5.54. For each shell, we analyze the spectra following 
the method reported in Sect. 3.1. Then we calculated the outflow 
properties as a function of radius from the launching region. We 
compute the mass outflow rate for each shell following Lutz et al. 
(2020) as 


Vout,she M. out,shell 


xR (8) 


Moutshell = 


where AR = 609 pc, being the size of each shell, Vout,shen is the 
velocity of the outflow computed as in Eq. 3, and Moutshen is 
the mass in the shell calculated as using Eq. 1. The outflow Ha 
luminosity was corrected for dust obscuration using the Balmer 
decrement computed in each shell. 

In this case, the mass outflow rate represents the amount of 
gas passing through each projected shell. We show the mass out- 
flow rate map in Fig. 14, in which each shell is color-coded by 
the mass outflow rate. 

In Fig. 15 we show the outflow velocity (computed using Eq. 
3), the mass outflow rate, the momentum rate, and the kinetic rate 
for each radial shell. We report all the quantities, assuming either 
no obscuration (Ay = 0, blue points) or the inferred obscuration 
based on the Balmer decrement in each shell (Ay + 0, purple 
points). In the first two rings and the last one, the ratio between 
the broad Ha and H8 components is lower than the theoretical 
ratio of 2.86, so we assumed Ay = 0. 

The observed outflow velocity is higher in the two rings clos- 
est to the launching region, then it decreases and stays almost 
constant up to 4 kpc. This can also be seen in the outflow maps 
from the spaxel-by-spaxel fitting in Fig. 6 and B.1. The mass 
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Fig. 15. Outflow velocity (top left), mass outflow rate (top right), momentum rate (bottom left), and kinetic rate (bottom right) in each radial 
shell as a function of radius from the central component (see 14). Blue and purple points indicate whether the outflow mass was calculated by 
assuming Ay = 0 or the Ay from the outflow Balmer decrement in each ring, respectively. The half-blue, half-purple points have measured Ay = 0. 
The purple and blue shaded regions in top-right panel show the 1o region around the mean outflow rate computed with Ay + 0 and Ay = 0, 


respectively. 


outflow rate for each shell instead reaches its maximum between 
1.4 and 3.5 kpc from the launching region. 

The momentum and kinetic rate radial profiles show a steady 
increase with radius up to 3.5 kpc in both quantities in the case 
in which we corrected for the dust attenuation of the outflow. 

If we assume that each shell has moved at a constant velocity 
Vout in the plane of the sky since its launch, then the gas in each 
shell represents the outflow ejected at a time T = Vout/l, where l 
is the intrinsic, not projected distance from the launching point 
to where it is observed now. With this assumption, we deduce 
that Mout was higher in the past since the mass outflow rate is 
higher at larger radii. 

Assuming that the outflow has been launched at a time tT we 
can compare the outflow observed at a radius R, with the SFR in 
the central clump at the time T = Vout Xx cos(i)/R, where Vout is the 
de-projected outflow velocity, R is the distance along the plane 
of the sky and i is the assumed inclination angle of the outflow 
with respect to the plane perpendicular to the line of sight. We 
performed a SED fitting using Bagpipes with the same param- 
eters as reported in Sect. 5, but with finer SFH bins in the last 20 
Myr to properly characterize the variation of the SFR in the last 
few Myrs when the outflow was launched. In Fig. 16 we plot the 


mass loading factor 7 as a function of lookback time, defined as 
j= Mion ou(R) /SFR(t), where t and R are related by the equa- 
tion reported above. For Vout we assumed 450 km sl, which is 
the average outflow velocity across the different shells. Varying 
the inclination angle i between 30 and 70 degrees we can observe 
that the outflows must have been stronger in the past, reaching 
values of 7 z 1 for i > 50. This analysis confirms what already 
observed by Herrera-Camus et al. (2021) with the [Cu] neutral 
outflow, the mass loading factor of the outflow is low consider- 
ing the entire galaxy (Nion ~ 3%), but becomes significant when 
considering only the central component. 


8. Discussion 


Based on previous observations, HZ4 was considered a typical 
high-z turbulent rotating disk. Our new JWST/NIRSpec IFS data 
instead reveal a clumpy structure, a different morphology de- 
pending on the tracer UV continuum emission, [On]45007 À and 
other rest-frame optical emission lines, [Cu] or dust continuum, 
with each clump having different SFHs and different properties, 
and an asymmetric velocity field. These new observations hint 
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Fig. 16. Mass loading factor of the ionized gas defined as Mionou/ SFR. 
as function of time for a range of assumed inclination angles. 


at this galaxy being a merger between at least two galaxies with 
different velocities, instead of a rotating disk. 

The presence of additional broader components, statistically 
needed to reproduce the observed line profile, also implies that 
there are non-circular motions both in the rest-frame optical 
lines, which trace the ionized gas, and in the re-analyzed [Cn] 
emission line, tracing mainly the neutral gas and confirming 
the previous results by Herrera-Camus et al. (2021). Following 
Herrera-Camus et al. (2021) we interpret these broad compo- 
nents as mainly due to outflows, but we discuss below different 
scenarios that can create the observed broad features (see Fig. 
17). 


— Different outflows in different directions. 
This scenario is supported by the observation of a star form- 
ing and obscured clump, which thanks to the SF feedback 
is able to launch a blue-shifted, ionized gas outflow with 
the energetics that we observe. On the other hand, the red- 
shifted neutral outflow can be launched by the same region 
or other star-forming clumps (like the [Oni] bright north- 
ern clump), triggered by the merger event. The ionized and 
molecular outflows are almost co-spatial, but the [Cr] out- 
flow extends in the northwest direction, while the ionized 
outflow extends in the northeast direction from what we as- 
sumed as the launching region. However, the Xsrg is high 
enough to allow the launch of these outflows from the same 
regions. We note that the observation of the redshifted out- 
flow only in [Cr] does not exclude the possibility of it also 
hosting ionized gas. If the receding side of the outflow lies 
behind the disk with respect to our line of sight, that side of 
the ionized outflow could be obscured by the dust in the sys- 
tem in the optical emission lines observed by NIRSpec, but 
still be observed in [Cu] with ALMA. In that case, we would 
only be able to detect the blueshifted approaching side. 
— Galactic fountain. 

The redshifted [Cu] emission can also be explained as a 
galactic fountain. The ionized wind is not fast enough to 
escape the potential well of the galaxy as shown in Sect. 
7.3, hence eventually the expelled material may cool down 
and then fall back onto the galaxy itself creating the galactic 
fountain (Leroy et al. 2015). It is expected that the majority 
of the gas launched outwards by SF-driven winds eventually 
will fall back as galaxy fountains (Bregman 1980; Oppen- 
heimer et al. 2010; Brook et al. 2012; Übler et al. 2014). 
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Given the low outflow velocity observed in comparison to 
the escape velocity, eventually this will happen, but it is un- 
clear whether we are witnessing it now in these observations. 
Unfortunately, the presence of a galactic fountain is hard to 
prove even in the local Universe and requires a detailed kine- 
matic model of the galaxy and the outflow to prove or dis- 
prove its presence (e.g., Yuan et al. 2023). 

- [Cu] inflow or merger-induced flows and [Om] outflow. 
The accretion of cold gas onto galaxies is one of the main 
mechanisms of galaxy growth and evolution. Inflow signa- 
tures are hard to observe and study at high redshift, and are 
usually seen thanks to absorption features (Rubin et al. 2012; 
Bouché et al. 2013; Herrera-Camus et al. 2020), or thanks to 
hints of radial motion via kinematic analysis (Arribas et al. 
2023; Genzel et al. 2023; Übler et al. 2024). 

If the most probable scenario is that the [Om] emission traces 
outflows, the [Cri] emission instead can be associated with 
the inflow of cold atomic and molecular gas from the CGM 
and IGM. As shown in Fig. 11, the orientation and location 
of the [Cu] broad component is aligned with the filamentary 
structure seen in the narrow component flux maps, both in 
the natural and in the Briggs cubes, which can resemble an 
accreting filament. 

With this scenario the accreting gas must be already enriched 
in [Cu], and thus cannot be pristine gas, hence the gas could 
originate from the galaxy, going back to the fountain sce- 
nario. 

Although streams of inflowing gas may be present in the 
system, it is unlikely that they dominate the observed broad 
line profiles, which imply motions of ~ 300 km s~!. While 
these velocities are typical of SF-driven outflows, they are 
too large for the inflow velocities that are expected to be a 
fraction of the virial velocity (vj; = VGMpw/ry, ~ 180 
km s^! ; Goerdt & Ceverino 2015). 

Another possibility is that the broad wings observed may be 
explained as due to gravitational interaction (1.e., tidal tails) 
between these galaxies in the process of merging (e.g., Baron 
et al. 2024). 


The present NIRSpec and ALMA observations describe 
a complex system formed by three merging galaxies, where 
likely gravitational interactions, inflowing gas, and outflows are 
present, with the latter dominating the large velocities implied 
by the broad component of the emission line profiles. A better 
knowledge of its actual geometry would allow us to further con- 
strain the properties of the system and analyze the relative role 
of these processes. 


9. Conclusions 


In this work, we have presented the JWST NIRSpec/IFS obser- 
vations of the galaxy HZ4 at z = 5.544. This galaxy is the highest 
redshift galaxy for which we have evidence of a neutral outflow 
traced by broad [Cu] emission lines (Herrera-Camus et al. 2021). 
We have analyzed both the low- and high-spectral resolution dat- 
acubes from JWST NIRSpec IFS to study for the first time the 
rest-frame UV and optical spectrum of this source. Moreover, we 
combined our results with a re-analysis of the [Cu] ALMA data 
to constrain the multiphase outflow properties of this galaxy. 
In particular, our main results are: 


— The JWST high spatial resolution observations reveal that 
the galaxy is not a regular rotating disk as inferred from 
ALMA lower resolution observations, but rather it is formed 
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Fig. 17. Cartoon illustrating three different scenarios for the explanation of the observations of the broad lines. In all panels, we draw in red the 
[Cu] non-circular components, and in blue the ionized gas emitting the broad rest-frame optical lines. In the left panel, we represent the multi-phase 
outflow scenario in which the ionized broad line is mainly emitted by the approaching side of the outflow launched by the SF clump. The outflow 
traced by the [Cu] broad component is emitted by outflowing gas on the receding side. In the central panel, we show the galactic fountain scenario, 
in which the ionized gas is still tracing the outflow, but the [Cir] emission is tracing the gas that cannot escape from the potential well of the galaxy, 
then it cools down and falls back onto the galaxy. The right panel instead depicts the scenario in which the redshifted [Cr] emission is tracing the 
inflow of cold gas from the CGM or IGM onto the galaxy, or is tracing gas stripped by the outer region of the northern and southern galaxies due 


to gravitational interaction. 


by three components. The north and the south components 
have different properties, such as stellar ages and star forma- 
tion histories. 

— We observe a spatially resolved ionized outflow, that, to- 
gether with the observation of the outflow in [Cn], allows 
us to probe the SF-driven multi-phase outflow for the first 
time in the early Universe. The ionized outflow extends over 
— 4 kpc from the launching site. 

— We constructed the ionized outflow history by analyzing the 
outflow properties in radial shells from the launching site. 
We observe that the velocity of the outflow decreases with 
increasing radius, while the mass outflow rate increases. We 
also computed the mass loading factor history of the outflow 
and found that the outflow was stronger in the past, reaching 
values of 7 ~ 1 — 2 at 5-20 Myrs, varying with the outflow 
inclination. 

— We find that in this z ~ 5.5 source the mass of the SF-driven 
outflow is dominated by the colder gas phases, in agreement 
with other low-z studies and a handful of studies up to z ~ 3. 


The synergy between ALMA and JWST observations has al- 
lowed us to probe the multi-phase ISM and outflow properties 
up to z ~ 5.5, but more observations and statistics are needed 
to probe the role of outflows in shaping galaxy properties in the 
early Universe. Increasing the number of high-redshift galaxies 
observed with multi-wavelength studies that target different out- 
flow phases (ionized, neutral atomic, molecular) will help us un- 
derstand the composition of outflows and constrain the role of 
outflows in galaxy evolution. 
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Appendix A: Other sources in the FOV 


While analyzing the R100 cube we detected another source in 
the field of view, which is distant ~0.8” from the northern bright 
clump of HZA. The source is located at RA = 09h58m28.55s, 
Dec = 02d03m06.69s, and it is undetected in HST or ALMA. 
The continuum emission of the source is negligible but is visi- 
ble only when integrating between 1.77 and 1.84 um as shown 
in the top panel of Fig. A.1. In the bottom panel of Fig. A.1 we 
show the spectrum of the source extracted from a circular aper- 
ture of radius 0.2" (4 spaxels) highlighted as the pink circle in the 
map. The spectra show two bright emission lines, that we iden- 
tify as the [Om]445007,4959À, HB complex and the Ha and the 
[N1]446548,6584À complex, which are unresolved in the prism 
spectrum. We performed a fit of the emission lines by using the 
same procedure as reported in 3.1 and implying a single Gaus- 
sian component for each emission line. This places the galaxy at 
a redshift of z = 2.607 + 0.003. 

We also estimate the galaxy properties by fitting the inte- 
grated spectrum with Bagpipes. 

The best-fit Bagpipes spectrum is shown in the bottom 
panel of Fig. A.1 as the red line. The galaxy has a stellar mass of 
log (M, /Mo) = 7.2 x 0.1. 
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Fig. A.1. Top panel: R100 cube integrated between 1.77 and 1.84 um. 
Bottom panel: In black we plot the spectrum extracted from the area 
highlighted as the magenta circle in the panel above. In red we show the 
best-fitting result from the Bagpipes fitting. 


Appendix B: Outflow velocity 


In this Sect. we show the outflow velocity map computed as Eq. 
3 for each spaxel and we compare it with the Xs rr map (see also 
Fig. 8). 


Appendix C: Outflow history 


In this Sect. we show the individual fit of each spectrum ex- 
tracted from each ring for the analysis of the outflow history in 
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Fig. B.1. Outflow velocity map defined as vou = |AVnarow broad] + 27 out 
(see also Sect. 7). We show as black dashed contours the Uspr map 
contours (see also Fig. 8). 


Fig. C.1. Ring 1 represents the closer one to the launching site, 
Ring 7 is the farther. We show the [Om], H6 complex. 
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Fig. C.1. Zoom in on the [Om] and H6 complex of each ring. Red and blue dashed lines are the narrow and outflow best-fit models. Please note 
that the displayed range on the y axis is the same for all the panels. 
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